####################################
# HTE effect RDD
####################################

rm(list=ls())

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)
library(gridExtra)

################
# Prepare data 
################

# read data
d = read.dta13("~/Dropbox/Anti-american attitudes/01_data/clean_data/trump_election_data_clean_final.dta")   
names(d)

# bandwidths
h_us_trust_bin = rdbwselect(d$us_trust_bin,d$score)$bws[1]
h_us_trust_bin

# covariates
table(d$female)
table(d$education)
table(d$age)
table(d$left_president)

# high school or more
d$highschool = 0
d$highschool[d$education>3]=1
table(d$highschool)

# older 50
d$more50 = 0
d$more50[d$age>49]=1
table(d$more50)

# interactions
d$int_treat_fem = d$treatment*d$female
d$int_treat_high = d$treatment*d$highschool
d$int_treat_more50 = d$treatment*d$more50
d$int_treat_leftpres = d$treatment*d$left_president
d$int_treat_south = d$treatment*d$south

############################
# US trust binary male
############################

# female
table(d$female)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$female + d$int_treat_fem + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$female + d$int_treat_fem + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$female + d$int_treat_fem + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_male = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") + 
  ggtitle("Male") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

#########################################
# US trust binary difference female male
#########################################

# female
table(d$female)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$female + d$int_treat_fem + d$score + as.factor(d$pais), weights = kweights)))[4, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$female + d$int_treat_fem + d$score + as.factor(d$pais), weights = kweights)))[4, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$female + d$int_treat_fem + d$score + as.factor(d$pais), weights = kweights)))[4, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_diff_female_male = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Difference female male") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

############################
# US trust binary female
############################

# male
d$male = 1 - d$female
d$int_treat_male = d$treatment*d$male

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$male + d$int_treat_male + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$male + d$int_treat_male + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$male + d$int_treat_male + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_female = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Female") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

#################################
# US trust binary no highschool
#################################

# highschool
table(d$highschool)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$highschool + d$int_treat_high + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$highschool + d$int_treat_high + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$highschool + d$int_treat_high + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_nohigh = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("No high-school") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

######################################################
# US trust binary difference highschool nohighschool
######################################################

# highschool
table(d$highschool)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$highschool + d$int_treat_high + d$score + as.factor(d$pais), weights = kweights)))[4, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$highschool + d$int_treat_high + d$score + as.factor(d$pais), weights = kweights)))[4, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$highschool + d$int_treat_high + d$score + as.factor(d$pais), weights = kweights)))[4, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_diff_high_nohigh = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Difference high-school no-high-school") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

################################
# US trust binary highschool
#################################

# no highschool
d$nohighschool = 1 - d$highschool
d$int_treat_nohigh = d$treatment*d$nohighschool
table(d$nohighschool)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$nohighschool + d$int_treat_nohigh + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$nohighschool + d$int_treat_nohigh + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$nohighschool + d$int_treat_nohigh + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_high = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("High school") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

#####################################
# US trust binary less than 50 years
#####################################

# More than 50
table(d$more50)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$more50 + d$int_treat_more50 + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$more50 + d$int_treat_more50 + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$more50 + d$int_treat_more50 + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_less50 = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Less than 50 years old") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

#########################################################
# US trust binary difference more than and less than 50
#########################################################

# More 50
table(d$more50)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$more50 + d$int_treat_more50 + d$score + as.factor(d$pais), weights = kweights)))[4, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$more50 + d$int_treat_more50 + d$score + as.factor(d$pais), weights = kweights)))[4, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$more50 + d$int_treat_more50 + d$score + as.factor(d$pais), weights = kweights)))[4, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_diff_more50_less50 = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Difference more and less than 50 years old") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

################################
# US trust binary more than 50
#################################

# less than 50
d$less50 = 1 - d$more50
d$int_treat_less50 = d$treatment*d$less50
table(d$int_treat_less50)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$less50 + d$int_treat_less50 + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$less50 + d$int_treat_less50 + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$less50 + d$int_treat_less50 + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_more50 = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("More than 50 years old") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

############################
# US trust binary non-left
############################

# left president
table(d$left_president)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$left_president + d$int_treat_leftpres + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$left_president + d$int_treat_leftpres + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$left_president + d$int_treat_leftpres + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_nonleft = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") + 
  ggtitle("Non-left-wing president") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

###########################################
# US trust binary difference left non-left
###########################################

# female
table(d$left_president)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$left_president + d$int_treat_leftpres + d$score + as.factor(d$pais), weights = kweights)))[4, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$left_president + d$int_treat_leftpres + d$score + as.factor(d$pais), weights = kweights)))[4, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$left_president + d$int_treat_leftpres + d$score + as.factor(d$pais), weights = kweights)))[4, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_diff_left_nonleft = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Difference left non-left presidents") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

############################
# US trust binary left
############################

# male
d$nonleft_president = 1 - d$left_president
d$int_treat_nonleftpres = d$treatment*d$nonleft_president

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$nonleft_president + d$int_treat_nonleftpres + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$nonleft_president + d$int_treat_nonleftpres + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$nonleft_president + d$int_treat_nonleftpres + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_left = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Left-wing presidents") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

############################
# US trust binary central
############################

# female
table(d$south)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$south + d$int_treat_south + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$south + d$int_treat_south + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$south + d$int_treat_south + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_central = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") + 
  ggtitle("Central America and the Caribbean") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

###############################################
# US trust binary difference south and central
################################################

# central
table(d$south)

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$south + d$int_treat_south + d$score + as.factor(d$pais), weights = kweights)))[4, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$south + d$int_treat_south + d$score + as.factor(d$pais), weights = kweights)))[4, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$south + d$int_treat_south + d$score + as.factor(d$pais), weights = kweights)))[4, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_diff_south_central = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("Difference South and Central America and the Caribbean") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

############################
# US trust binary south
############################

# central
d$central = 1 - d$south
d$int_treat_central = d$treatment*d$central

# loop
kweights = 0
pe = rep(NA,32)
se = rep(NA,32)
pv = rep(NA,32)
n  = rep(NA,32)

for (i in 1:32) {
  
  kweights = kernelwts(d$score, 0, bw = (i+5), kernel = "triangular")
  pe[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$central + d$int_treat_central + d$score + as.factor(d$pais), weights = kweights)))[2, 1]
  se[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$central + d$int_treat_central + d$score + as.factor(d$pais), weights = kweights)))[2, 2]
  pv[i] = coef(summary(lm(d$us_trust_bin ~ d$treatment + d$central + d$int_treat_central + d$score + as.factor(d$pais), weights = kweights)))[2, 4]
  n[i] = sum(kweights!=0 & !is.na(d$us_trust_bin))
}

days = c(6:37)
ts = pe/se
lower95 = pe - se*1.960
upper95 = pe + se*1.960
lower90 = pe - se*1.650
upper90 = pe + se*1.650
color = 0
results = data.frame(days,pe,se,ts,pv,n,lower90,upper90,lower95,upper95,color)
results$color[(results$days==round(h_us_trust_bin))==T] = 1
results$color[results$color==1]="red"
results$color[results$color==0]="black"
table(results$color)
head(results)

f_south = ggplot(results, aes(x=days, y=pe)) + geom_pointrange(aes(ymin=lower90,ymax=upper90,color=color),size=1.5,position = position_dodge(0.5)) + geom_pointrange(aes(ymin=lower95,ymax=upper95,color=color)) + scale_x_continuous(breaks=c(5,10,15,20,25,30,35,40)) + scale_y_continuous(breaks=c(-0.2,-0.1,0,0.1,0.2)) + coord_cartesian(ylim=c(-0.2,0.2)) + geom_hline(yintercept=0) + xlab("Bandwidth (absolute values)") + ylab("RDD estimates") + scale_colour_manual(values = c("black" = "black", "red" = "red")) + theme(axis.text=element_text(size=20)) + theme_grey(base_size = 20) + theme(legend.position="none") +
  ggtitle("South America") + theme(plot.title = element_text(hjust = 0.5)) + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank()) + theme(plot.title = element_text(size = 20, face = "bold"))

################################
# Plots
#################################

grid.arrange(f_female, f_male, f_diff_female_male, ncol = 1)

grid.arrange(f_high, f_nohigh, f_diff_high_nohigh, ncol = 1)

grid.arrange(f_more50, f_less50, f_diff_more50_less50, ncol = 1)

grid.arrange(f_left, f_nonleft, f_diff_left_nonleft, ncol = 1)

grid.arrange(f_south, f_central,f_diff_south_central,  ncol = 1)


